# Code for Graphs and Tables West European Politics paper - UK  analysis R and R

library(tidyverse)
#library(Hmisc)
library(lubridate)
library(broom)
library(gganimate)
library(transformr)
library(gtools)
library(jtools)
library(huxtable)
library(estimatr)
library(texreg)
library(modelsummary)


# Import datasets ---------------------------------------------------------

# Note new dataset is all local authorities and up to December - but limiting to existing time periods to be
# consistent with other analysis


uk_final <- read_csv("uk_final_dec20.csv", guess_max =  20000)

uk_final <- uk_final %>% 
  filter(date<as.Date("2020-08-08"))

weekend_df <- data.frame(date = unique(uk_final$date)) %>%
  mutate(day = weekdays(date)) %>%
  mutate(weekend= case_when(
    day %in% c("Saturday", "Sunday") ~ 1,
    date==as.Date("2020-04-10") ~ 1,
    date==as.Date("2020-04-13") ~ 1,
    date==as.Date("2020-05-08") ~ 1,
    date==as.Date("2020-05-25") ~ 1,
    date==as.Date("2020-08-31") ~ 1,
    TRUE ~ 0))

# UK Dragon Plot -------------------------------------------------------------

uk_final %>% 
  ggplot(aes(x = date, y = workplace, group=location, alpha = ALL_AGES))+
  geom_jitter(aes(color=as.factor(final_remain>50)))+
  scale_color_manual(values=c('blue','yellow'))+
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
 # geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
 #geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
  #geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  ylab("Change in Workplace Activity from Baseline")+xlab("Date")+
  annotate(geom="text", x=as.Date("2020-03-07"), y=12, label="Pre-Lockdown", color="red", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-03-14"), y=18, label="Full Lockdown", color="red", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-04-27"), y=10, label="Recreation Loosened", color="blue", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-06-03"), y=10, label="Retail Opened", color="blue", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-06-24"), y=18, label="'Independence'\n Day", color="blue", size = 2.5)+
 # annotate(geom="text", x=as.Date("2020-09-01"), y=10, label="Rule of Six", color="red", size = 2.5)+
  #annotate(geom="text", x=as.Date("2020-10-18"), y=12, label="Second Lockdown", color="red", size = 2.5)+
 # annotate(geom="text", x=as.Date("2020-11-22"), y=8, label=" End of Lockdown", color="blue", size = 2.5)+
  theme_classic()+
  theme(legend.position = "none")


ggsave("uk_dragon_randr.pdf", width = 8, height = 5)
ggsave("uk_dragon_randr.jpg", width = 8, height = 5)



# UK workplace regression plots -------------------------------------------



# Begin with ISCO one

coefs_map_workplace <- map_df(unique(uk_final$date), function(date_i) {
  filtered_data <- filter(uk_final, date == date_i) %>%
    lm(data =., workplace ~ final_remain+isco_ratio+final_gdp_cap+log_av_price+pop_under_30+pop_over_70+log(density)+I(region)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_workplace <- coefs_map_workplace %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("final_remain", "isco_ratio", "final_gdp_cap", "log_av_price", "pop_over_70", "log(density)")) %>% 
  mutate(term = case_when(term=="final_remain" ~ "Remain Vote",
                          term=="isco_ratio" ~ "Occupational Ratio",
                          term=="final_gdp_cap" ~ "GDP per Capita",
                          term=="log_av_price" ~ "Log House Prices",
                          term == "pop_over_70" ~ "Population > 70",
                          term == "log(density)" ~ "Log Density"))


coefs_map_workplace %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dotted", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
  geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  xlab("Date")+ylab("Day by Day Coefficients (DV: Workplace Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +labs(title="Workplace")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("uk_coefs_isco.pdf", width = 8, height = 5)
ggsave("uk_coefs_isco.png", width = 8, height = 5)

# WITH SECTORS ADDED AND ISCO REMOVED

coefs_map_workplace_2 <- map_df(unique(uk_final$date), function(date_i) {
  filtered_data <- filter(uk_final, date == date_i) %>%
    lm(data =., workplace ~ final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+pop_under_30+pop_over_70+log(density)+I(region)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})


coefs_map_workplace_2 <- coefs_map_workplace_2 %>%
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>%
  filter(term %in% c("final_remain", "prim_sec_sector", "pop_over_70", "log_av_price", "high_skill_serv_sector", "log(density)")) %>%
  mutate(term = case_when(term=="final_remain" ~ "Remain Vote",
                          term=="prim_sec_sector" ~ "Primary/ Sec Sectors",
                          term == "pop_over_70" ~ "Population > 70",
                          term=="log_av_price" ~ "Log House Prices",
                          term == "high_skill_serv_sector" ~ "High Skilled Services",
                          term == "log(density)" ~ "Log Density"))

coefs_map_workplace_2 %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dotted", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
  geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  xlab("Date")+ylab("Day by Day Coefficients (DV: Workplace Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +labs(title="Workplace")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("uk_coefs_sectors.pdf", width = 8, height = 5)
ggsave("uk_coefs_sectors.png", width = 8, height = 5)

# WITH COVID ADDED AND ISCO AND SECTORS

coefs_map_workplace_3 <- map_df(unique(uk_final$date), function(date_i) {
  filtered_data <- filter(uk_final, date == date_i) %>%
    lm(data =., workplace ~ final_remain+week_cases_rate+isco_ratio+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+pop_under_30+pop_over_70+log(density)+I(region)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})


coefs_map_workplace_3 <- coefs_map_workplace_3 %>%
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>%
  filter(term %in% c("final_remain", "week_cases_rate", "final_gdp_cap", "log_av_price", "pop_over_70", "log(density)")) %>%
  filter(!(term == "week_cases_rate" & date < as.Date("2020-03-11"))) %>% 
  mutate(term = case_when(term=="final_remain" ~ "Remain Vote",
                          term=="week_cases_rate" ~ "Weekly Cases Rate",
                          term=="final_gdp_cap" ~ "GDP per Capita",
                          term=="log_av_price" ~ "Log House Prices",
                          term == "pop_over_70" ~ "Population > 70",
                          term == "log(density)" ~ "Log Density"))


coefs_map_workplace_3 %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dotted", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
  geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  xlab("Date")+ylab("Day by Day Coefficients (DV: Workplace Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +labs(title="Workplace")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("uk_coefs_cases.pdf", width = 8, height = 5)
ggsave("uk_coefs_cases.png", width = 8, height = 5)


# UK residential regression plots -----------------------------------------



coefs_map_residential <- map_df(unique(uk_final$date), function(date_i) {
  filtered_data <- filter(uk_final, date == date_i) %>%
    lm(data =., residential ~ final_remain+isco_ratio+final_gdp_cap+log_av_price+pop_under_30+pop_over_70+log(density)+I(region)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_residential <- coefs_map_residential %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("final_remain", "isco_ratio", "final_gdp_cap", "log_av_price", "pop_over_70", "log(density)")) %>% 
  mutate(term = case_when(term=="final_remain" ~ "Remain Vote",
                          term=="isco_ratio" ~ "Occupational Ratio",
                          term=="final_gdp_cap" ~ "GDP per Capita",
                          term=="log_av_price" ~ "Log House Prices",
                          term == "pop_over_70" ~ "Population > 70",
                          term == "log(density)" ~ "Log Density"))


coefs_map_residential %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dotted", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
  geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  xlab("Date")+ylab("Day by Day Coefficients (DV: Residential Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +labs(title="Residential")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("uk_coefs_res_isco.pdf", width = 8, height = 5)
ggsave("uk_coefs_res_isco.png", width = 8, height = 5)

## With sectors

coefs_map_residential_2 <- map_df(unique(uk_final$date), function(date_i) {
  filtered_data <- filter(uk_final, date == date_i) %>%
    lm(data =., residential ~ final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+pop_under_30+pop_over_70+log(density)+I(region)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_residential_2 <- coefs_map_residential_2 %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("final_remain", "prim_sec_sector", "pop_over_70", "log_av_price", "high_skill_serv_sector", "log(density)")) %>%
  mutate(term = case_when(term=="final_remain" ~ "Remain Vote",
                          term=="prim_sec_sector" ~ "Primary/ Sec Sectors",
                          term == "pop_over_70" ~ "Population > 70",
                          term=="log_av_price" ~ "Log House Prices",
                          term == "high_skill_serv_sector" ~ "High Skilled Services",
                          term == "log(density)" ~ "Log Density"))


coefs_map_residential_2 %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dotted", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
  geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  xlab("Date")+ylab("Day by Day Coefficients (DV: Residential Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +labs(title="Residential")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("uk_coefs_res_sector.pdf", width = 8, height = 5)
ggsave("uk_coefs_res_sector.png", width = 8, height = 5)

# WITH COVID ADDED AND ISCO AND SECTORS

coefs_map_residential_3 <- map_df(unique(uk_final$date), function(date_i) {
  filtered_data <- filter(uk_final, date == date_i) %>%
    lm(data =., residential ~ final_remain+week_cases_rate+isco_ratio+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+pop_under_30+pop_over_70+log(density)+I(region)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})


coefs_map_residential_3 <- coefs_map_residential_3 %>%
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>%
  filter(term %in% c("final_remain", "week_cases_rate", "final_gdp_cap", "log_av_price", "pop_over_70", "log(density)")) %>%
  filter(!(term == "week_cases_rate" & date < as.Date("2020-03-11"))) %>% 
  mutate(term = case_when(term=="final_remain" ~ "Remain Vote",
                          term=="week_cases_rate" ~ "Weekly Cases Rate",
                          term=="final_gdp_cap" ~ "GDP per Capita",
                          term=="log_av_price" ~ "Log House Prices",
                          term == "pop_over_70" ~ "Population > 70",
                          term == "log(density)" ~ "Log Density"))


coefs_map_residential_3 %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dotted", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
  geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  xlab("Date")+ylab("Day by Day Coefficients (DV: Residential Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +labs(title="Residential")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("uk_coefs_res_cases.pdf", width = 8, height = 5)
ggsave("uk_coefs_res_cases.png", width = 8, height = 5)


# UK regressions ----------------------------------------------------------


# UK regression analysis

uk_final <- uk_final %>% 
  mutate(after_march_16 = as.numeric(date>as.Date("2020-03-16")),
         after_march_22 = as.numeric(date>as.Date("2020-03-22")),
         after_may_12 = as.numeric(date>as.Date("2020-05-12")),
         after_june_14 = as.numeric(date>as.Date("2020-06-14")),
         after_july_3 = as.numeric(date>as.Date("2020-07-03"))
       #  after_sep_13 = as.numeric(date>as.Date("2020-09-13")),
      #   after_nov_4 = as.numeric(date>as.Date("2020-11-04")),
       #  after_dec_1 = as.numeric(date>as.Date("2020-12-01"))
      ) %>% 
  group_by(location) %>% 
  mutate(lag_workplace = lag(workplace),
         l2_workplace = lag(workplace, 2),
         l3_workplace = lag(workplace, 3),
         l4_workplace = lag(workplace, 4),
         l5_workplace = lag(workplace, 5),
         l6_workplace = lag(workplace, 6)) %>% 
  ungroup()

library(estimatr)
s3<-lm_robust(data = uk_final, workplace ~ lag_workplace+  final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                final_remain*after_june_14+final_remain*after_july_3+
                #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = region)

s3_l2<-lm_robust(data = uk_final, workplace ~ lag_workplace+ l2_workplace+ final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                   final_remain*after_june_14+final_remain*after_july_3+
                   #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                   final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density),  fixed_effects = region)

s3_l6<-lm_robust(data = uk_final, workplace ~ lag_workplace+ l2_workplace+ l3_workplace + l4_workplace + l5_workplace + l6_workplace+ final_remain*after_march_16+
                   final_remain*after_march_22+final_remain*after_may_12+final_remain*after_june_14+final_remain*after_july_3+
                   #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                   final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density),  fixed_effects = region)



s3_1<-lm_robust(data = uk_final, workplace ~ lag_workplace+ final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                  final_remain*after_june_14+final_remain*after_july_3+ 
                  #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                  final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)


s3_1_l2<-lm_robust(data = uk_final, workplace ~ lag_workplace+l2_workplace+ final_remain*after_march_16+final_remain*after_march_22+
                     final_remain*after_may_12+final_remain*after_june_14+final_remain*after_july_3+
                     #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                     final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)

s3_1_l6<-lm_robust(data = uk_final, workplace ~ lag_workplace+l2_workplace+l3_workplace + l4_workplace + l5_workplace + l6_workplace+
                     final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+final_remain*after_june_14+
                     final_remain*after_july_3+
                     #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                     final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)


r1 <- lm_robust(data = uk_final, workplace ~ lag_workplace+
                  (final_remain)*after_march_16+(final_remain)*after_march_22+
                  (final_remain)*after_may_12+(final_remain)*after_june_14+
                  (final_remain)*after_july_3,
                  #(final_remain)*after_sep_13+(final_remain)*after_nov_4+(final_remain)*after_dec_1 
                  fixed_effects = location)


r2 <- lm_robust(data = uk_final, workplace ~ lag_workplace+
                  (final_remain)*after_march_16+(final_remain)*after_march_22+
                  (final_remain)*after_may_12+(final_remain)*after_june_14+
                  (final_remain)*after_july_3
                  #(final_remain)*after_sep_13+(final_remain)*after_nov_4+
                  #(final_remain)*after_dec_1
                , fixed_effects = date+location)

r1_full <- lm_robust(data = uk_final, workplace ~ lag_workplace+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_16+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_22+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_may_12+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_june_14+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_july_3
                      # (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_sep_13+
                      #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_nov_4+
                      #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_dec_1
                      ,                      fixed_effects = location)


r2_full <- lm_robust(data = uk_final, workplace ~ lag_workplace+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_16+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_22+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_may_12+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_june_14+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_july_3
                      #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_sep_13+
                      #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_nov_4+
                      #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_dec_1
                      ,                   fixed_effects = date+location)


coefs_s3 <- list( "lag_workplace" = "Lag workplace" , 
                  "final_remain" = "Remain vote" ,  "isco_ratio" = "ISCO Ratio"  , "final_gdp_cap" = "GDP per cap","log_av_price" =  "Log (House prices)"  ,
                  "pop_under_30" = "Pop < 30" ,  "pop_over_70" = "Pop > 70" , "log(density)" = "Log (Density)" ,  "after_march_16" = "After Mar 16",
                  "after_march_22" = "After Mar 22" , "after_may_12" = "After May 12" ,  "after_june_14" = "After June 14", "after_july_3" = "After July 3" ,
                  "after_sep_13" = "After Sep 13", "after_nov_4" = "After Nov 4", "after_dec_1" = "After Dec 1",
                  "final_remain:after_march_16" = "Remain * After Mar 16", "final_remain:after_march_22" = "Remain * After Mar 22",
                  "final_remain:after_may_12"  = "Remain * After May 12", "final_remain:after_june_14" = "Remain * After June 14" ,
                  "final_remain:after_july_3" =  "Remain * After July 3" , "final_remain:after_sep_13" =  "Remain * After Sep 13",
                  "final_remain:after_nov_4" =  "Remain * After Nov 4", "final_remain:after_dec_1" = "Remain * After Dec 1" )

library(texreg)
# This table has Model 1 


texreg(list(s3, s3_1, r1, r2, r1_full, r2_full), include.ci = FALSE, omit.coef = "(l2)|(l3)|(l4)|(l5)|(l6)", custom.coef.map = coefs_s3,
       custom.gof.rows = list( "Region Effects"= c("Y", "Y", "N", "N", "N", "N"), "Date Effects" = c( "N", "Y", "N","Y", "N", "Y"),
                               "Locality Effects" = c("N", "N", "Y", "Y", "Y", "Y"), "All Interactions" = c("N", "N", "N", "N", "Y", "Y"))
, file="UK_workplace_RandR.tex", caption = "Workplace Activity in the UK", 
caption.above = TRUE, include.rsquared=FALSE, include.rmse=FALSE, label = "table_uk", booktabs = TRUE
)


uk_final_no_weekend <- uk_final %>% 
  filter(weekend==0) 



s3w<-lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+  final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                final_remain*after_june_14+final_remain*after_july_3+
                 #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                 prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = region)

s3_l2w<-lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+ l2_workplace+ final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                   final_remain*after_june_14+final_remain*after_july_3+
                    #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                    prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density),  fixed_effects = region)

s3_l6w<-lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+ l2_workplace+ l3_workplace + l4_workplace + l5_workplace + l6_workplace+ final_remain*after_march_16+
                   final_remain*after_march_22+final_remain*after_may_12+final_remain*after_june_14+
                    #final_remain*after_july_3+
                   #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                    prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density),  fixed_effects = region)



s3_1w<-lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+ final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                  final_remain*after_june_14+final_remain*after_july_3+ 
                   #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                   prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)


s3_1_l2w<-lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+l2_workplace+ final_remain*after_march_16+final_remain*after_march_22+
                     final_remain*after_may_12+final_remain*after_june_14+final_remain*after_july_3+
                     #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                      prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)

s3_1_l6w<-lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+l2_workplace+l3_workplace + l4_workplace + l5_workplace + l6_workplace+
                     final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+final_remain*after_june_14+
                     final_remain*after_july_3+
                      #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                      prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)


r1w <- lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+
                  (final_remain)*after_march_16+(final_remain)*after_march_22+
                  (final_remain)*after_may_12+(final_remain)*after_june_14+
                  (final_remain)*after_july_3
                   #(final_remain)*after_sep_13+(final_remain)*after_nov_4+
                  #(final_remain)*after_dec_1
                 , fixed_effects = location)


r2w <- lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+
                  (final_remain)*after_march_16+(final_remain)*after_march_22+
                  (final_remain)*after_may_12+(final_remain)*after_june_14+
                  (final_remain)*after_july_3
                 #+(final_remain)*after_sep_13+(final_remain)*after_nov_4+
                  #(final_remain)*after_dec_1
                 , fixed_effects = date+location)

r1w_full <- lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_16+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_22+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_may_12+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_june_14+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_july_3
                       # (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_sep_13+
                        #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_nov_4+
                        #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_dec_1
                        ,                      fixed_effects = location)


r2w_full <- lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+
                   (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_16+
                   (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_22+
                   (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_may_12+
                   (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_june_14+
                   (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_july_3
                   #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_sep_13+
                   #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_nov_4+
                   #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_dec_1
                     ,                   fixed_effects = date+location)

coefs_s3 <- list( "lag_workplace" = "Lag workplace" , 
                  "final_remain" = "Remain vote" , "prim_sec_sector" = "Manual Employment", "high_skill_serv_sector" = "High Skill Services" , "isco_ratio" = "ISCO Ratio"  , "final_gdp_cap" = "GDP per cap","log_av_price" =  "Log (House prices)"  ,
                  "pop_under_30" = "Pop < 30" ,  "pop_over_70" = "Pop > 70" , "log(density)" = "Log (Density)" ,  "after_march_16" = "After Mar 16",
                  "after_march_22" = "After Mar 22" , "after_may_12" = "After May 12" ,  "after_june_14" = "After June 14", "after_july_3" = "After July 3" ,
                  "after_sep_13" = "After Sep 13", "after_nov_4" = "After Nov 4", "after_dec_1" = "After Dec 1",
                  "final_remain:after_march_16" = "Remain * After Mar 16", "final_remain:after_march_22" = "Remain * After Mar 22",
                  "final_remain:after_may_12"  = "Remain * After May 12", "final_remain:after_june_14" = "Remain * After June 14" ,
                  "final_remain:after_july_3" =  "Remain * After July 3" , "final_remain:after_sep_13" =  "Remain * After Sep 13",
                  "final_remain:after_nov_4" =  "Remain * After Nov 4", "final_remain:after_dec_1" = "Remain * After Dec 1" )


# This table has Model 1 
texreg(list(s3w, s3_1w, r1w, r2w, r1w_full, r2w_full), include.ci = FALSE, omit.coef = "(l2)|(l3)|(l4)|(l5)|(l6)", custom.coef.map = coefs_s3,
       custom.gof.rows = list( "Region Effects"= c("Y", "Y", "N", "N", "N", "N"), "Date Effects" = c( "N", "Y", "N","Y", "N", "Y"),
                              "Locality Effects" = c("N", "N", "Y", "Y", "Y", "Y"), "All Interactions" = c("N", "N", "N", "N", "Y", "Y"))
       , file="UK_workplace_weekdays_RandR.tex", caption = "Workplace Activity in the UK (No Weekends)", 
       caption.above = TRUE, include.rsquared=FALSE, include.rmse=FALSE, label = "table_uk", booktabs = TRUE
)

# Predicted dragon



## Appendix Materials

### UK Retail and Recreation graphs

coefs_map_retail <- map_df(unique(uk_final$date), function(date_i) {
  filtered_data <- filter(uk_final, date == date_i) %>%
    lm(data =., retail_rec ~ final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+pop_under_30+pop_over_70+log(density)+I(region)) %>%
    broom::tidy() %>%
    rename(coef = estimate, se = std.error) %>%
    mutate(date = date_i)
})

coefs_map_retail <- coefs_map_retail %>% 
  mutate(upper = coef+2*se,
         lower = coef-2*se) %>% 
  filter(term %in% c("final_remain", "prim_sec_sector", "pop_over_70", "log_av_price", "high_skill_serv_sector", "log(density)")) %>%
  mutate(term = case_when(term=="final_remain" ~ "Remain Vote",
                          term=="prim_sec_sector" ~ "Primary/ Sec Sectors",
                          term == "pop_over_70" ~ "Population > 70",
                          term=="log_av_price" ~ "Log House Prices",
                          term == "high_skill_serv_sector" ~ "High Skilled Services",
                          term == "log(density)" ~ "Log Density"))


coefs_map_retail %>% 
  left_join(weekend_df) %>%
  filter(weekend==0)  %>%
  ggplot(aes(x = date, y = coef))+
  geom_line()+
  geom_hline(yintercept = 0, linetype = "dotted", colour = "black") +
  geom_ribbon(aes(ymin=lower, ymax = upper, alpha=0.3), color="lightgray")+
  #lockdown
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
 # geom_vline(xintercept =as.Date("2020-09-14"), linetype="dashed" , color = "red")+
 # geom_vline(xintercept =as.Date("2020-11-05"), linetype="dashed", color = "red")+
 # geom_vline(xintercept =as.Date("2020-12-02"), linetype="dashed", color = "blue")+
  xlab("Date")+ylab("Day by Day Coefficients (DV: Retail and Recreation Activity)")+
  theme_minimal()+
  theme(legend.position = "none") +labs(title="Retail and Recreation")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  facet_wrap(~term, scales = "free_y")

ggsave("uk_coefs_retail_sector.pdf", width = 8, height = 5)
ggsave("uk_coefs_retail_sector.png", width = 8, height = 5)



## UK residential regression table

uk_final <- uk_final %>% 
  mutate(after_march_16 = as.numeric(date>as.Date("2020-03-16")),
         after_march_22 = as.numeric(date>as.Date("2020-03-22")),
         after_may_12 = as.numeric(date>as.Date("2020-05-12")),
         after_june_14 = as.numeric(date>as.Date("2020-06-14")),
         after_july_3 = as.numeric(date>as.Date("2020-07-03")),
        # after_sep_13 = as.numeric(date>as.Date("2020-09-13")),
         #after_nov_4 = as.numeric(date>as.Date("2020-11-04")),
         #after_dec_1 = as.numeric(date>as.Date("2020-12-01"))
        ) %>% 
  group_by(location) %>% 
  mutate(lag_residential = lag(residential),
         l2_residential = lag(residential, 2),
         l3_residential = lag(residential, 3),
         l4_residential = lag(residential, 4),
         l5_residential = lag(residential, 5),
         l6_residential = lag(residential, 6)) %>% 
  ungroup()

uk_final_no_weekend <- uk_final %>% 
  filter(weekend==0)

s3w<-lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+  final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                 final_remain*after_june_14+final_remain*after_july_3+
               #+final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                 prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = region)

s3_l2w<-lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+ l2_residential+ final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                    final_remain*after_june_14+final_remain*after_july_3+
                  #+final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                    prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density),  fixed_effects = region)

s3_l6w<-lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+ l2_residential+ l3_residential + l4_residential + l5_residential + l6_residential+ final_remain*after_march_16+
                    final_remain*after_march_22+final_remain*after_may_12+final_remain*after_june_14+final_remain*after_july_3+
                    #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                    prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density),  fixed_effects = region)



s3_1w<-lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+ final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                   final_remain*after_june_14+final_remain*after_july_3+
                   #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                   prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)


s3_1_l2w<-lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+l2_residential+ final_remain*after_march_16+final_remain*after_march_22+
                      final_remain*after_may_12+final_remain*after_june_14+final_remain*after_july_3+
                      #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                      prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)

s3_1_l6w<-lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+l2_residential+l3_residential + l4_residential + l5_residential + l6_residential+
                      final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+final_remain*after_june_14+
                      final_remain*after_july_3+
                      #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                      prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density), fixed_effects = ~date+region)


r1w <- lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+
                   (final_remain)*after_march_16+(final_remain)*after_march_22+
                   (final_remain)*after_may_12+(final_remain)*after_june_14+
                   (final_remain)*after_july_3
                   #(final_remain)*after_sep_13+(final_remain)*after_nov_4+
                   #(final_remain)*after_dec_1
                   , fixed_effects = location)


r2w <- lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+
                   (final_remain)*after_march_16+(final_remain)*after_march_22+
                   (final_remain)*after_may_12+(final_remain)*after_june_14+
                   (final_remain)*after_july_3
                 #+(final_remain)*after_sep_13+(final_remain)*after_nov_4+
                   #(final_remain)*after_dec_1
                 , fixed_effects = date+location)

r1w_full <- lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_16+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_22+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_may_12+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_june_14+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_july_3
                        #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_sep_13+
                       # (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_nov_4+
                       #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_dec_1,
                      , fixed_effects = location)


r2w_full <- lm_robust(data = uk_final_no_weekend, residential ~ lag_residential+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_16+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_march_22+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_may_12+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_june_14+
                        (final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_july_3
                        #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_sep_13+
                       #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_nov_4+
                        #(final_remain+prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density))*after_dec_1,
                      , fixed_effects = date+location)

coefs_s3 <- list( "lag_residential" = "Lag Residential" , 
                  "final_remain" = "Remain vote" , "prim_sec_sector" = "Manual Employment", "high_skill_serv_sector" = "High Skill Services" , "isco_ratio" = "ISCO Ratio"  , "final_gdp_cap" = "GDP per cap","log_av_price" =  "Log (House prices)"  ,
                  "pop_under_30" = "Pop < 30" ,  "pop_over_70" = "Pop > 70" , "log(density)" = "Log (Density)" ,  "after_march_16" = "After Mar 16",
                  "after_march_22" = "After Mar 22" , "after_may_12" = "After May 12" ,  "after_june_14" = "After June 14", "after_july_3" = "After July 3" ,
                  "after_sep_13" = "After Sep 13", "after_nov_4" = "After Nov 4", "after_dec_1" = "After Dec 1",
                  "final_remain:after_march_16" = "Remain * After Mar 16", "final_remain:after_march_22" = "Remain * After Mar 22",
                  "final_remain:after_may_12"  = "Remain * After May 12", "final_remain:after_june_14" = "Remain * After June 14" ,
                  "final_remain:after_july_3" =  "Remain * After July 3" , "final_remain:after_sep_13" =  "Remain * After Sep 13",
                  "final_remain:after_nov_4" =  "Remain * After Nov 4", "final_remain:after_dec_1" = "Remain * After Dec 1" )


# This table has Model 1 
texreg(list(s3w, s3_1w, r1w,  r2w, r1w_full, r2w_full), include.ci = FALSE, omit.coef = "(l2)|(l3)|(l4)|(l5)|(l6)", custom.coef.map = coefs_s3,
       custom.gof.rows = list( "Region Effects"= c("Y", "Y", "N",  "N", "N", "N"), "Date Effects" = c( "N", "Y", "N", "Y", "N", "Y"),
                               "Locality Effects" = c("N", "N", "Y", "Y",  "Y", "Y"), "All Interactions" = c("N", "N", "N", "N", "Y", "Y"))
       , file="UK_residential_weekdays_RandR.tex", caption = "Residential Activity in the UK", 
       caption.above = TRUE, include.rsquared=FALSE, include.rmse=FALSE, label = "table_res_uk", booktabs = FALSE
)


# Interaction with income - split sample by median GDP per capita

uk_final <- uk_final %>% 
  mutate(after_march_16 = as.numeric(date>as.Date("2020-03-16")),
         after_march_22 = as.numeric(date>as.Date("2020-03-22")),
         after_may_12 = as.numeric(date>as.Date("2020-05-12")),
         after_june_14 = as.numeric(date>as.Date("2020-06-14")),
         after_july_3 = as.numeric(date>as.Date("2020-07-03"))
         #  after_sep_13 = as.numeric(date>as.Date("2020-09-13")),
         #   after_nov_4 = as.numeric(date>as.Date("2020-11-04")),
         #  after_dec_1 = as.numeric(date>as.Date("2020-12-01"))
  ) %>% 
  group_by(location) %>% 
  mutate(lag_workplace = lag(workplace),
         l2_workplace = lag(workplace, 2),
         l3_workplace = lag(workplace, 3),
         l4_workplace = lag(workplace, 4),
         l5_workplace = lag(workplace, 5),
         l6_workplace = lag(workplace, 6)) %>% 
  ungroup()

uk_final_no_weekend <- uk_final %>% 
  filter(weekend==0) 

s3h<-r2w <- lm_robust(data = uk_final_no_weekend %>%  filter(final_gdp_cap>median(final_gdp_cap)), workplace ~ lag_workplace+
                        (final_remain)*after_march_16+(final_remain)*after_march_22+
                        (final_remain)*after_may_12+(final_remain)*after_june_14+
                        (final_remain)*after_july_3
                      #+(final_remain)*after_sep_13+(final_remain)*after_nov_4+
                      #(final_remain)*after_dec_1
                      , fixed_effects = date+location)

s3l<-r2w <- lm_robust(data = uk_final_no_weekend %>%  filter(final_gdp_cap<=median(final_gdp_cap)), workplace ~ lag_workplace+
                        (final_remain)*after_march_16+(final_remain)*after_march_22+
                        (final_remain)*after_may_12+(final_remain)*after_june_14+
                        (final_remain)*after_july_3
                      #+(final_remain)*after_sep_13+(final_remain)*after_nov_4+
                      #(final_remain)*after_dec_1
                      , fixed_effects = date+location)

summary(s3h)
summary(s3l)

coefs_s3 <- list( "lag_workplace" = "Lag Workplace" , 
                  "final_remain:after_march_16" = "Remain * After Mar 16", "final_remain:after_march_22" = "Remain * After Mar 22",
                  "final_remain:after_may_12"  = "Remain * After May 12", "final_remain:after_june_14" = "Remain * After June 14" ,
                  "final_remain:after_july_3" =  "Remain * After July 3" , "final_remain:after_sep_13" =  "Remain * After Sep 13",
                  "final_remain:after_nov_4" =  "Remain * After Nov 4", "final_remain:after_dec_1" = "Remain * After Dec 1" )

texreg(list(s3h, s3l), include.ci = FALSE, omit.coef = "(l2)|(l3)|(l4)|(l5)|(l6)", custom.coef.map = coefs_s3,
       custom.gof.rows = list( "Income Level"  = c("High", "Low"), "Date Effects" = c("Y", "Y"),
                               "Locality Effects" = c( "Y", "Y"))
       , file="UK_workplace_split_RandR.tex", caption = "Workplace Activity in the UK Split by Income", 
       caption.above = TRUE, include.rsquared=FALSE, include.rmse=FALSE, label = "table_res_uk", booktabs = FALSE
)


# Simulation plot for UK



s3w<-lm_robust(data = uk_final_no_weekend, workplace ~ lag_workplace+  final_remain*after_march_16+final_remain*after_march_22+final_remain*after_may_12+
                 final_remain*after_june_14+final_remain*after_july_3+
                 #final_remain*after_sep_13+final_remain*after_nov_4+final_remain*after_dec_1+
                 prim_sec_sector+high_skill_serv_sector+final_gdp_cap+log_av_price+isco_ratio+pop_under_30+pop_over_70+log(density)+ I(region))

fitted = predict(s3w, uk_final %>% filter(region!="Northern Ireland"))

uk_final2 <- cbind(uk_final, .fitted)

uk_final2 %>% 
  filter(weekend==0) %>% 
  ggplot(aes(x = date, y = .fitted, group=location, alpha = ALL_AGES))+
  geom_line(aes(color = as.factor(final_remain>50), group=location))+
  scale_color_manual(values=c('blue','yellow'))+
  geom_vline(xintercept =as.Date("2020-03-17"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-03-23"), linetype="dashed", color = "red")+
  geom_vline(xintercept =as.Date("2020-05-13"), linetype="dashed", color = "blue")+
  geom_vline(xintercept =as.Date("2020-06-15"), linetype="dashed" , color = "blue")+
  geom_vline(xintercept =as.Date("2020-07-04"), linetype="dashed" , color = "blue")+
  ylab("Predicted Change in Workplace Activity from Baseline")+xlab("Date")+
  annotate(geom="text", x=as.Date("2020-03-09"), y=10, label="Pre-Lockdown", color="red", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-03-17"), y=15, label="Lockdown", color="red", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-04-29"), y=10, label="Recreation Loosened", color="blue", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-06-05"), y=10, label="Retail Opened", color="blue", size = 2.5)+
  annotate(geom="text", x=as.Date("2020-06-25"), y=8, label="'Independence'\n Day", color="blue", size = 2.5)+
  theme_classic()+
  theme(legend.position = "none")

ggsave("simulation.pdf", width = 8, height = 5)
ggsave("simulation.pdf", width = 8, height = 5)
